﻿/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;

namespace dnAnalytics.Statistics.Distributions
{
    /// <summary>
    /// The Dirichlet distribution is a distribution over probability vectors. It is parameterized by a vector of
    /// non-negative numbers.
    /// </summary>
    /// <remarks>The distribution will use the <see cref="System.Random"/> by default. 
    /// Users can set the random number generator by using the <see cref="RandomNumberGenerator"/> property.
    /// 
    /// The statistics classes will check all the incoming parameters whether they are in the allowed
    /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
    /// to false, all parameter checks can be turned off.</remarks>
    public class Dirichlet
    {
        // The Dirichlet distribution parameters.
        private readonly double[] mAlpha;

        /// <summary>
        /// Constructs a new Dirichlet distribution.
        /// </summary>
        /// <param name="alpha">An array with the Dirichlet parameters.</param>
        public Dirichlet(double[] alpha)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(alpha);
            }

            mAlpha = alpha;
        }

        /// <summary>
        /// Constructs a new symmetric Dirichlet distribution.
        /// </summary>
        /// <param name="alpha">The value of each parameter of the Dirichlet distribution.</param>
        /// <param name="k">The dimension of the Dirichlet distribution.</param>
        public Dirichlet(double alpha, int k)
        {
            if (Control.CheckDistributionParameters)
            {
                if (alpha < 0.0)
                {
                    throw new ArgumentOutOfRangeException("alpha", "Alpha entries should be positive.");
                }

                if (k <= 0)
                {
                    throw new ArgumentOutOfRangeException("k", "The dimension of the Dirichlet distribution should be positive.");
                }
            }

            mAlpha = new double[k];
            for (int i = 0; i < k; i++)
            {
                mAlpha[i] = alpha;
            }
        }

        /// <summary>
        /// A string representation of the distribution.
        /// </summary>
        public override string ToString()
        {
            return "Dirichlet(Dimension = " + this.Dimension + ")";
        }

        /// <summary>
        /// The dimension of the Dirichlet distribution.
        /// </summary>
        public int Dimension
        {
            get { return mAlpha.Length; }
        }

        /// <summary>
        /// The parameters of the Dirichlet distribution.
        /// </summary>
        public double[] Alpha
        {
            get { return mAlpha; }
        }

        /// <summary>
        /// The sum of all alpha parameters.
        /// </summary>
        public double AlphaSum
        {
            get
            {
                double s = 0.0;
                for (int i = 0; i < mAlpha.Length; i++)
                {
                    s += mAlpha[i];
                }
                return s;
            }
        }

        /// <summary>
        /// Gets the mean.
        /// </summary>
        /// <value>The mean.</value>
        public double[] Mean
        {
            get { return mAlpha; }
        }

        /// <summary>
        /// Gets the variance.
        /// </summary>
        /// <value>The variance.</value>
        public double[] Variance
        {
            get
            {
                double s = this.AlphaSum;
                double[] v = new double[mAlpha.Length];
                for (int i = 0; i < mAlpha.Length; i++)
                {
                    v[i] = mAlpha[i]*(s - mAlpha[i])/(s*s*(s + 1.0));
                }
                return v;
            }
        }

        /// <summary>
        /// Gets or sets the random number generator.
        /// </summary>
        /// <value>The random number generator used to generate a random sample.</value>
        public System.Random RandomNumberGenerator { get; set; }

        /// <summary>
        /// Check the parameters of the distribution.
        /// </summary>
        /// <exception cref="ArgumentOutOfRangeException">If any of the components in the Dirichlet parameter vector are negative.</exception>
        public static void CheckParameters(double[] alpha)
        {
            for (int i = 0; i < alpha.Length; i++)
            {
                if (alpha[i] < 0.0)
                {
                    throw new ArgumentOutOfRangeException("alpha", "Alpha entries should be positive.");
                }
            }
        }

        /// <summary>
        /// Samples a Dirichlet distributed random vector.
        /// </summary>
        public double[] Sample()
        {
            return Sample(RandomNumberGenerator, mAlpha);
        }

        /// <summary>
        /// Samples a Dirichlet distributed random vector.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="alpha">The Dirichlet distribution parameter.</param>
        public static double[] Sample(System.Random rnd, double[] alpha)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(alpha);
            }

            int n = alpha.Length;
            double[] gv = new double[n];
            double sum = 0.0;
            for (int i = 0; i < n; i++)
            {
                gv[i] = Gamma.Sample(rnd, alpha[i], 1.0);
                sum += gv[i];
            }

            for (int i = 0; i < n; i++)
            {
                gv[i] /= sum;
            }

            return gv;
        }
    }
}